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ABSTRACT 

—An  implicit  form  of  upstream  differencing  is  developed  which 
is  decoupled  in  space  and  is  unconditionally  stable.  It  has  the 
same  numerical  diffusion  as  the  usual  explicit  form.  The  scheme 
is  especially  attractive  for  use  in  radiation  boundary  conditions 
for  semi-implicit  models. 
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This  note  examines  an  alternate  form  of  upstream  differencing  which  is 
implicit.  Although  it  could  be  used  as  a  scheme  for  advection  in  a  finite 
difference  model,  our  interest  is  motivated  by  a  need  for  an  implicit  radiation 
boundary  condition.  We  wish  to  use  a  Sommerfeld  radiation  condition  of  the  form 


u.  +  (u  +  c)  u  =0 
t  x 


where  u  is  a  velocity  component,  c  is  a  gravity  wave  speed,  and  the  subscripts  are 
time  and  space  derivatives  (see  Klemp  and  Lilly,  1978).  In  an  explicit  model  this 
can  be  implemented  by  employing  the  usual  explicit  form  of  upstream  differencing. 

In  a  model  where  gravity  waves  are  treated  Implicitly  it  is  likely  that  the  time 

step  would  violate  the  criterion  for  linear  stability  of  an  explicit  upstream 

scheme.  An  implicit  upstream  scheme  which  was  unconditionally  stable  would  be  of 
great  value  in  this  situation.  However,  the  model  design  requires  further  restric¬ 
tions  on  the  formulation  of  such  a  scheme,  if  it  is  to  be  useful  in  the  model. 
Namely,  it  must  be  possible  to  calculate  the  boundary  values  implicitly  prior  to  and 

independent  of  the  calculation  of  the  interior  values  at  the  new  time  level  (see 

Hurlburt  and  Thompson,  1980).  That  is,  the  implicit  scheme  must  be  decoupled  in 
space.  We  will  attempt  to  design  an  implicit  upstream  differencing  formula  for 
Equation  (1)  which  satisfies  all  the  preceding  requirements. 


To  examine  upstream  differencing  consider  the  advect ion-diffusion  equation 


ut  +  uux  -  Auxx  -  0 


which  we  difference  in  space  using  centered  schemes  which  are  second-order  accurate: 


ut  "  "  2 Ax 


(uj+i  •  vi)+  r?  (vi  -  2uj +  vi) 


l 


wnere  j  is  a  grid  index  and  Ax  is  a  grid  increment  for  u.  If  we  let  the  coefficient 
of  the  diffusion  term  equal  the  coefficient  of  the  advective  term,  then 


A  ■  %UjAx 


and  Equation  (3)  becomes 


ut  51  ■  £ 


(4) 


(5) 


Equation  (5)  is  the  usual  explicit  form  of  upstream  differencing  for  advection,  if  a 
forward  (Euler)  time  difference  is  used.  It  is  first -order  accurate  in  space 
and  time.  Equations  (3)  and  (4)  show  that  numerical  diffusion  is  implicit  in 
Equation  (5).  Table  1  shows  the  equivalent  eddy  viscosity  for  different  values  of  u 
and  Ax  using  (4). 


Table  1:  Numerical  eddy  viscosity  for  upstream  differencing 


A  u  Ax 

10®  cm^/s  100  cm/s  20  km 
10?  cm^/s  10  cm/s  20  km 
10®  cm^/s  1  cm/s  20  km 


In  mesoscale  ocean  eddy  problems  with  Ax  -  20  km  and  u^x  -100  cm/s,  typically 
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A  -  3x10  cm  8  has  been  used.  Hence,  upstream  differencing  is  highly  diffusive. 
It  is  less  diffusive  than  A  -  3x10®  cm^/s  when  Ax*20  km  only  where  |u|  <3  cm/s. 
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For  a  leapfrog  time  difference, the  advective  term  in  Equation  (3)  la  conditionally 
stable;  but  the  diffusive  term  is  unstable,  and  the  upstream  difference  is  unstable. 
For  an  Euler  time  difference,  the  advective  term  is  unstable  and  the  diffusion  term 
is  conditionally  stable.  The  upstream  difference  is  conditionally  stable  (At<Ax/u) 
for  u  >  0.  Note  that  forming  the  upstream  difference  from  Equation  (3)  requires 
that  Uj+i  and  Uj_j  be  at  the  same  time  level  in  both  the  advective  and  diffusive 
terms.  However,  this  is  not  required  for  uj.  Nor  does  it  require  uj  at  the  same 
time  level  as  uj+j  and  Uj_j.  Hence,  It  is  possible  to  use  leapfrog  for  the 
advective  term  and  Dufort  Frankel  for  the  diffusive  term.  In  this  case  (3)  becomes 


n+1 


n-1 


2AtUj 

IaT 


(  n  _n  \  .  2AtA  f  n 

Ivi '  ui-ij +  v* 


-  u. 


n+1 


j+i  j 


-  U. 


n-1 


Vi) 


(6) 


where  n  is  a  time  step  index  and  At  is  a  time  increment.  If,  as  before,  we  use 
Equation  (4),  then 


n+1 


=  u, 


n-1 


2u^At 

Ax 


(7) 


Equation  (7)  is  an  implicit  upstream  differencing  formula  which  is  decoupled  in 
space  as  desired  for  the  radiation  boundary  condition  in  our  semi-implicit  ocean 
model.  Although  second  order  accurate  in  time,  Equation  (7)  is  first  order  in  space 
and  has  the  same  numerical  diffusion  as  Equation  (5).  Also,  Equation  (7)  is  a  three 
time  level  scheme  whereas  Equation  (5)  required  only  two  time  levels.  Thus, 

Equation  (7)  has  a  real  practical  advantage  over  Equation  (5)  only  if  it  is 
unconditionally  stable.  Since  in  Equation  (6)  the  advective  term  is  conditionally 
stable  and  the  diffusive  term  is  unconditionally  stable,  we  anticipate  that  Equation 
(7)  is  at  least  conditionally  stable.  Since  Equation  (7)  is  implicit  and  since  the 
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stability  of  the  diffusive  tern  in  Equation  (3)  was  dominant  when  u  >  0,  we  have 
some  hope  for  unconditional  stability  for  Equation  (7).  We  will  test  this  hypothe¬ 
sis  using  a  linear  stability  analysis. 

Consider  the  liner ized  advection  equation 


u 


t 


-  -Uu 


x 


(8) 


represented  by  the  difference  equation.  Equation  (7).  We  linearize  Equation  (7)  by 
letting  U*"Uj  in  the  coefficient  of  the  space  difference.  Also  let 


n 


u 


J 


ikjAx 

u  e 
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(9) 


and  substitute  in  Equation  (7)  to  get 


n+1 


u 


n-1 


UAt 

Ax 


n+l 


+  u 


n-1 


-  2u  e 
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-ikAx 


(10) 


Let 


B  -  UAt/Ax 


(ID 


the  Courant  number,  and  let 


a  **  kAx 


(12) 
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Then  Equation  (10)  can  be  written 


2B  -la  ,  (  1-B  \ 

Un+1  =  1+B  Une  +  KTti)  Vl 


We  need  to  write  Equation  (13)  so  that  the  coefficients  form  an  amplification 
matrix,  G,  l.e.. 


u  , ,  =  Gu 
n+1  n 


Thus,  Equation  (13)  becomes 


n+1  1+B 


2B  _lor  1-B 


1+B  n 


If  the  spectral  radius,  0  of  G  <  1,  then  un  is  bounded  as  n-*-<»,  where  o  -  max 
(Am)  and  the  Xm  are  the  eigenvalues  of  the  matrix.  We  calculate  the  eigen¬ 
values  using  the  characteristic  equation  obtained  from  det  j XI— G | —0 ,  where  1  Is  the 
identity  matrix.  The  characteristic  equation  for  Equation  (15)  is 


.2  2B  -ia  .  1-B 

A  '  TTb*  x  "  TTb 


which  has  roots 


_  -ia  f  / 
B  e  ±  ( 

1+B  V 


-ia  \2 


+  1± 
J  1+B 
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The  value  of  IM  is  maximum  when  e“ia  -  +1.  The  use  of  both  is  redundant,  so  we 
consider  only  e”*a  ■  1.  For  the  positive  root 


tnAx  |  A  |  ■  1 


(18) 


for  the  negative 


max  |  \ 


(19) 


From  Equations  (11),  (18),  and  (19),  Equation  (7)  is  unconditionally  stable  for 
U  >  0  (upstream  differencing)  and  is  unstable  for  U  <  0  (downstream  differencing). 

Thus,  an  implicit  upstream  difference  has  been  formulated  which  is 
unconditionally  stable  and  is  decoupled  in  space.  It  meets  all  the  requirements  for 
implementation  in  a  radiation  boundary  condition  in  our  semi-implicit  ocean  model. 
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